General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



K77-1J74b 


(•.«ASA-C*-l49 31b) SiNGUiiitilli CufiFUTAIiLNS 
(cdcnegie lost, or lech.) J7 y h C AoJ/MF 
AO 1 CSCL 12A 

Uaclas 

G3/64 57b7d 


SINGULARITY COMPUTATIONS 


J. L. Swedlow 


Report SM-76-6 


October 1976 


This paper was prepared as an invited lecture at the 
International Symposium on Innovative Numerical Methods 
in Applied Engineering Science 
Versailles, 23-27 May 1977 




Department of Mechanical Engineering 
Carnegie Institute of Technology 
Carnegie-Mellon University 
Pittsburgh, Pennsylvania 


A 

DEC 1976 



% 

received 

NASKStl fcClLKY ; 
input BRANCH 

c 


l 

5 



SINGULARITY COMPUTATIONS 


J. L. Swedlow 

Carnegie- Me lion University 

One of the intriguing - and sometimes perplexing - classes of problems in 
mechanics involves singularities in otherwise smooth fields. Examples abound 
from the Joukowsky airfoil to Kelvin’s problem and hardly need recounting 
here. Where an analytical solution is available, numerics may be sim- 
plified or avoided altogether. In other instances, numerical analysis is 
necessary, and properties of the singularity are then inferred from tab- 
ular information. These data are typically sparse or inaccurate in the 
immediate region of the singularity, or the numerical technique affords poor 
resolution, or some other impediment is encountered in establishing fully 
the result required. 

Needed is an approach inherently untroubled by such shortcomings, one 
that indicates the structure of a singularity directly. In particular, it 
would be useful to have both .the radial and angular (polar or spherical) 
distributions of the field quantities delineated as explicitly as is prac- 
ticable, together with some measure (s) of the intensity of the singularity. 

In this paper, we suggest such an approach, based on recent development of 
numerical methods for elasto-plastic flow. This approach is patently appli- 
cable to other problems in solid mechanics and, without much effort, lends 
itself to certain types of heat flow, fluid motion, and the like. 

Analytic solutions to classical problems in mechanics where a singularity 
occurs are divided, for our purposes, into two classes. In the first, one 
variable or set of variables is finite at the origin of the singularity, but 



2 


a second set - typically gradients - is not. Prime examples are displace- 
ments and strains at re-entrant corners, temperature and heat flux at an 
abrupt change in surface insulation. The second class involves singularities 
where none of the quantities is everywhere finite as, for example, in 
Kelvin’s problem. This type of solution has proven essential to development 
of numerical procedures such as the boundary integral method and it thereby 
deserves close attention; because concern here rests with finite element 
methods in solid mechanics, where boundedness of displacements is necessary, 
we limit attention to just the first class of problem. 


PLANAR ELASTICITY 

Singular behavior at re-entrant comers in classical planar elasticity 
(plane stress, plane strain) was fully articulated by Professor M. L. Williams [1] 
long before the widespread use of finite elements. His basic result in 1952 
can be interpreted as the sum of two series, one of which provides the singu- 
larity in strain and stress, if such behavior exists for the geometry and 
boundary conditions prescribed. The second series gives regular results 
which, in the finite element sense, provide the components of rigid motion, 

’’linear" displacements or "constant" strains*, and increasingly higher order 
terms. This latter component of Williams's solution may be shown to be 
equivalent to the interpolation functions used commonly in regular finite 
elements. 

There have been many developments reported in the literature to incor- 
porate the first part of Williams's solution into finite elements for the 
specific case of a crack; see, e.g., [2,3]. Professor P. M. Quinlan has 
also used these functions to enhance an edge-function form of analysis [4], 
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and Dr. B. Gross and his colleagues have published a variety of solutions 
using collocation methods [5]. In all these instances, and in others where 
extraction of a stress intensity is the objective [6], it is essential to 
have at the outset a certain foreknowledge of the structure of the singularity. 
That is, were it not known from Williams - or an equivalent source - what to 
expect, none of these procedures could have been implemented. 

Where information as to local behavior is desired, it is obvious that 
the details must be properly characterized either in advance or as an- integral 
part of the (numerical) analysis. The case of planar elasticity is well in 
hand and no further primary findings are to be anticipated. For the purpose 
of our discussion, however, we make the supposition that the structure of 
the singularity - if, indeed, one were to exist - is not known and then 
enquire what steps might be taken to reveal that information. Knowing the 
basis of Williams’s work as well as his results should then provide guidance 
for devising a numerical procedure. 

Let us presume that, in some specific problem of interest, displacements 
in the vicinity of a suspected singularity behave in the following manner. 

In addition to rigid motion and linear variations in the displacement field 
which produce the familiar ’’constant" strain's, the displacements along a 
ray from the origin (i.e., the point of singularity) behave as 

u 'v p^ * (1) 

where p is a (linear) radial coordinate, and q is an undetermined exponent. 

Let us furthermore focus on the case where 0 < q < 1 to provide the type of 
problem in which we are interested, i.e., singular gradient (s) of u. It is 
noted that no regular element will produce the response shown in (1); while 
for some range of p there may be a correspondence of regular element behavior 
and (1), the similarity is fortuitous and cannot be relied upon. 



Our task now becomes that of finding q in (1), pertinent to some specific 
problem. In fact, the approach employed here addresses this task and is con- 
ceptually no different than that used in finite element analysis - we simply 
go a step further. An element is defined and its displacement field pre- 
sumed as is usual, with terms of the form (1) included. Strains are com- 
puted, and the potential energy is defined and minimized with respect to 
nodal displacements; the result is the familiar stiffness equation. In a 
particular problem, however, solution strategy extends beyond solving the 
stiffness equation: potential energy is simultaneously minimized with 

respect to the exponent in (1) to achieve the full result. 

To be more specific, we assume the re-entrant comer to be surrounded 
by an array of sectors which together comprise a special element, as in 
Figure 1. Arbitrarily here, each sector has five nodes as sketched and we 
assume the cartesian components of displacement take the form 

u = u + A.x + C y + (E+F0)p^ cos 0 - (G+H0)p^sin 0 

1 (2a) 

v = v + D^x + B y + (E+F6)p^sin0 + (G+H6)p^cos 9 

« 

Accordingly, the polar components of displacement are 

u = u cos 0 + v sin 0 + ^-(A-.+B-Jp + ^-(A,-B,)p cos 20 
o o 2 1 1 2 1 1 

+ 7^ C i +D i^ p sin 20 + ( E+F0 )p q (2b) 

v. = -u sin 0 + v cos 0 - )p + i(C, +D.)p cos 20 

O O O L 1 1 l i I 

- y(A^-B^)p sin 20 + (G+H0)p^ 

Clearly, u q ,v q , and (D^-Cj) are rigid motions; further, u^v^A^.H are 
coefficients to be determined (in terms of nodal displacements). The process 
of relating these coefficients to the nodal values, although tedious, is 
straightforward and the result is easily tabulated [7]. 
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The "special" aspect of the element thus outlined is evident in the 
terms in (2). This simple representation is not more than an extension of 
the well-known constant strain element (CSE) and is thus but one of many 
possible configurations that one might contemplate for use in a given situ- 
ation. For the present purpose, we take the element in Figure 1 to be em- 
bedded in an array of "regular" or unmodified CSEs which make up the struc- 
ture to be analyzed as, for example, is shown in Figure 2. Some other 
possibilities are outlined in the Appendix. 

The result of having evaluated the coefficients in (2) in terms of 
nodal displacements may be written 

{u} = [o(p f e)]{u> (3) 

where {u} represents the two components in (2a) and {u} is a vector of 
nodal displacement components. The matrix [a] is the set of interpolation 
functions; note that q appears only in [a]. Following standard procedure, 
we next compute from (3) the strains 
fe. 


e y } = U) = [s(p, 9) I (u) " 

v 

xyj 


(4) 


as a precursor to computing the potential energy of the entire element 
assembly. In this context, note that (3) and (4) pertain to a typical 
sector of the special element, so that the potential energy will involve 
a summation of the contributions from each sector. 

Since 

la ) 


xy; 


{ o } = [M] {e} 



we have the potential energy for the typical sector as 


U 

s 




[B] T [M] [8]pdpde{u} 
[a] T {t }dS 


a 

where {t} is the traction vector specified on some part of the sector's 
boundary S. Summing terms of the form (5) over the sectors, we arrive at 
the special element's potential energy 


( 5 ) 



in which many of the traction terms obviously cancel. We next add the poten- 
tial energies of each of the regular elements U r to arrive at the total value 
for the system: 

U = U + l U (6) 

e r r 

To the extent that the surfaces of the re-entrant corner are traction-free 
or, alternatively, that their nodal displacement components are specified, 
the only uncancelled contributions to the overall system's traction are 
precisely where tractions are specified, normally far from the re-entrant 
comer. 

We are now prepared to minimize the functional U with respect to nodal 
displacements, arriving at the familiar stiffness equation. In the present 
development, the entries in the stiffness matrix take the same form for 
regular (CSE) elements as has been shown in many places. The contribution 
from each sector of the special element derives from the first integral in (5) 
and, together, a stiffness for the special element may readily be identified. 
Even with this relation, however, it is clear that no constraint on the expo- 



nent q has obtained; indeed, solving the stiffness equation could go forward 
with virtually any value for q. Were q somehow known, however, the computation 
at this stage is unremarkable except for (what appears to be) some negligible 
inconrpatabilities between the special and regular elements. 

The constraint on q, so far lacking, is obtained by minimizing the func- 
tional U with respect to the exponent itself. Were we to set 3U/3q = 0, how- 
ever, the result would be an hideously non-linear equation to be handled along 
with the stiffness equation. It is apparent from (6) that minimization of U 
with respect to q is equivalent to minimization of U e , since U r is not func- 
tionally dependent upon q. Hence, if the result of minimizing U with respect 
to all nodal displacements (u) is the familiar stiffness equation 
[K]{u> * {T} 

we have the simultaneous statement 

U g * minimum with respect to q (7b) 

Together, (7a) and (7b) pose the problem fully for prescribed nodal forces {T}. 

It is well at this point to comment on the connection between this formu- 
lation and that employed by Williams [1], so that the equivalence emerges . It 
will be recalled that Williams employed the Airy stress function x(p*8) in 
the form 

X(p.e) = p X+1 F(0;A) (8) 

and required x(p» 0 ) to be, a biharmonic function. This led to an ordinary 
differential equation for F (8 ; A) whose solution gave the proper dependence 
on 0, with four constants of integration, but left A unspecified. Williams 
argued that A > 0 is necessary to provide finite displacements as p + 0 
(since the displacement components were demonstrably 0(p X )), and that 
0 < A < 1 would lead to singularities in strain and stress. To determine 
A, Williams invoked sets of boundary conditions at 0 = 0 Q and 0 = -0 Q (see 



Figure 1); these sets of conditions were limited to the surfaces of the re- 
entrant comer being clamped (zero displacements) or free (zero tractions). 

The four conditions (two at 8 q two at -e ) for the four constants of integra- 
tion are thus homogeneous; the necessary condition for their solution yields 
an eigenequation for A. The result, especially for a crack (0 Q - it), has 
been used repeatedly in mechanics as noted above. 

If, alternatively, Williams had begun by writing displacements in the form 

u p = p f(6; A) 

v 0 - p x g(e; a) 

and required the stresses derived from (9) to satisfy the equations of equili- 
brium, the eventual result would be identical to his published findings. What 
is done here, of course, follows an alternate procedure which results in one 
significant difference. The representation (2) is not required to satisfy 
equilibrium at an arbitrary point ( p , 0) but is forced to do so over a finite 
region, here, the sector. That is, the result of minimizing U imposes 
equilibrium element- (or sector-) wise and not point-wise. 

As to the boundary conditons, there is no difference . Were (2) inserted 
into an analytic development ’based on minimum potential energy, the same boundary 
conditions as used by Williams would obtain. Thus the exponent q in (2) is 
subject to precisely the same constraint Williams obtained for A in (8), or 
would have obtained for A in (9). 

What the present formulation then provides is an approximate statement 
of interior equilibrium and an effectively exact boundary condition. The 
interior approximation is perforce tailored to the form of the assumed dis- 
placement components, whether they be (2) or (9) or forms suggested in the 
Appendix. The boundary constraint, however, is identical to that in Williams’s 
eigenequation. It may be noted further that Williams' s argument for the 
finiteness of displacement as p -*■ 0 is here replaced by boundedness of 



of the potential energy, a condition familiar to numerical analysts and, in 
certain respects, more easily treated. 

Returning to the problem statement (7), we comment on a solution tactic. 
Both (7a) and (7b) must be satisfied, and it appears straightforward to pro- 
ceed on an iterative basis. Let us assume a reasonable starting value for q and 
solve the linear equations in (7a) in the same manner as used in dealing with 
standard finite element problems. Arriving at a set of nodal displacements, 
the function is minimized with respect to the exponent q. Note that since 
nodal displacements are fixed, S a in (5) is null and the minimization pro- 
ceeds on the first integral in (5), summed over all sectors of the special 
element. A new value of q is found and (7a) is solved again. The process 
continues until both (7a) and (7b) are satisfied to whatever degree of pre- 
cision is appropriate to the computation and the computer involved. Having 
thus converged, the solution, in terms of (U } and q, is in hand. Straight- 
forward data reduction will provide both the full structure of the singularity, 
at least within the precision of (2), and the intensity of the singularity 
for the problem of interest. 

SOME SIMPLE EXTENSIONS 

The formulation outlined above may be extended to allied problems. The 
inclusion of body forces due, for example, to gravitational or centrifugal 
loading (as in a spinning disc) is effected merely through adding appropriate 
terms to the potential energy in (5) and (6) . Thermal excitation is included 
in an analogous manner. Dynamic behavior is modeled by replacing the theorem 
of minimum potential energy by Hamilton's principle. Such extensions in 
ordinary planar elasticity are theoretically well founded [8] and operational 
in a number of existing codes. 



Another type of extension to the foregoing development is equally ob- 
vious. Material anisotropy is incorporated via simple changes in the matrix 
[M]. Spatial variation in material properties may also be incorporated by 
appropriate alterations to [M], although the analyst should take note of the 
spatial gradients in [M] when designing an element map. With due caution, 
then, problems of the sort considered by Hein and Erdogan [9], for example, 
may be treated by finite elements. 


NON-LINEARITIES 

To this point, we have done little more that show an alternate technique 
for replicating the basic information contained in Williams’s eigenanalyses. 
While his findings give a basis for demonstrating and substantiating the 
present approach, its utility derives from circumstances where an eigen- 
analyses aoes not exist or is available solely through highly idealized 
modeling. 

Non-linear behavior is a case in point. If we consider first material 
non-linearity due to yield, we observe that the only analytic result available 
is the so-called HRR model [10]. This situation pertains to planar behavior, 
as does Williams's work, but is limited in certain respects. It admits plastic 
deformations only, it idealizes the material, and results presently available 
pertain only to a crack and not to the general re-entrant comer. Since in 
large measure the effect of yielding (or plastic flow, or non-linear material 
response) is to alter the initial singularity, it would be most useful to 
know how the change proceeds from the outset. That is, there may be signi- 
ficant technological interest in the process (as well as its rate of progress) 
whereby the material goes from the one limiting case described by Williams 
to the other limit characterized by the HRR model. 
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We have given considerable attention to this issue, and results are 
just now coming to hand. A preliminary discussion appears elsewhere [11], 
and further documentation is anticipated, e.g., [12]. For the prusent dis- 
cussion a brief outline of the formulation is in order. The steps described 
above are followed except that (3) is replaced by 
{6u}* [a(p,0)]{5u) 

where the 6 signifies an increment in each of the various displacement compo- 
nents. Then (4) becomes 
{6c} * [0(p,9}]{6u> 

and since the flow rule for elasto-plastic flow is written 
{60} - [M] { 6 e } 

one needs only an analogue to (5) to carry the analysis through. The re- 
quired theorem :s in fact available [l.*>] and we write for (7) 

[K]{6u} = {6T1 (10a) 

U ^ = minimum with respect to q (10b) 

The problem (10) is to be solved successively for the incremental values of 
the displacement components, not their accumulated or cur.-ent values. This 
problem is linear in the increment - in the sense that (7) is linear - and 
the procedure for its solution is established [14]. 

Certain matters relating to implementation are to be noted. One must 
choose the radial extent of the special element, the number of sectors used, 
and the refinement within each sector needed to effect proper quadrature. To 
investigate the~e matters, extensive evaluation of the code was performed using 
elementary solutions (as in [1]) as a basis. No clear criteria for sizing 
the element, i.e., fixing p . emerged except where a relationship could be 
established a priori by a given set of loading conditions. Circumferential 



behavior, on the other hand, could easily be seen to improve with increasing 
numbers of sectors, at least for the simple formulation (2). It does appear 
from work in progress [15] that a more refined representation behaves in the 
same manner, reaching acceptable performance with a modest number of sectors 
10]. For the element described in (2), however, we chose 48 sectors for 
-ir < 0 < n as the best trade-off between angular resolution and storage require- 
ments*. Quadrature was evaluated using a number of techniques; while the first 
integral in (5) is always bounded, it contains non-analytic function of p 
which impede a formal prediction of the behavior of various methods. In 
the event, we chose Gaussian quadrature with three angular positions in 
0 1 < ® < ®2 an< * seven i n 0 < p < 2p e . Good accuracy was obtained for 
various values of q in (?); other combinations may become preferred in other 
models, e.g., as described in the Appendix. 

The residual issue of sizing the radial extent of the special element 
was resolved empirically. For the center-cracked configuration of [11], 
we examined the elastic potential energy of the test specimen as •' was 
reduced, finding that a stationary value obtained for p e ^ 0 (crack length/ 100). 
While surely this result reflects both the modeling in (2) and numerical pre- 
cision of the computation itself, it also provides confidence in the element 
mapping. Having thus arrived at a suitable element array, an elasto-plastic 

analysis for a configuration reported jarlier [16] proceeded; in this manner 

# 

distinct results for the same problem were available to ensure that the special 

element computations could be corioborated. Other problems have since been 

considered and results are to be reported shortly [12,15]. It is of interest 

here, however, to note one or two aspects of the solution in [11]. Data 

♦Actually, analyses were performed for a symmetric configuration so that 24 
sectors were employed for the half specimen. 



were presented for four load levels: 

purely elastic response; 

yield detected just beyond the special element, 
denoted as load step 38; 

yield extended through the cross-section, 
denoted as load step 73; and 

average applied stress exceeded the yield point, 
denoted as load step 93. 

The radial variation of the octahedral or an effective stress, normalized 
on the respective value of yield, is shown in Figure 3 for one angular posi- 
tion . Radial variation of u - u and v, normalized on the uniform far-field 

o 

extension A, is shown in Figure 4*. Angular distributions of the same quan- 
tities appear in Figures 5 and 6. Note that the elastic results which exhibit 
high gradients are smoothly described; at high yield levels, roughness develops 
in some of the data. Nonetheless, one is able even with a crude representation 
of the type given in (2) to infer a fair sense of the structure of the crack 
tip’s singularity as yield proceeds. Incidentally, it may also be remarked 
that the analysis compares favorably to the experimentally observed specimen 
behavior for which a companion analysis was performed earlier [16]. 

As a second item, geometric non-linearity may be incorporated in the 
analysis. While the one analytic solution [17] to the problem confirms a 
localized singularity, it is necessarily confined to a specialized material 
representation. The obviotis issues then become, for a more arbitrary material 
characterization, the degree to which the sharp corner blunts and the size 
scale over which this event occurs as loading proceeds. Dr. J. R. Osias has 
investigated this matter using conventional elements in his original Eulerian 
formulation [18]. More recently, he has reformulated the problem for a special 
•Nomenclature is the same as that used in (2a); for reasons of symmetry v q = 0. 
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element of the type considered here, necessarily using a Lagrangian coor- 
dinate frame [19]. In this manner, the radial coordinate whose exponent is 
to be determined is readily identified. 

It is useful to touch briefly on the field problem. The initial domain 
B q is bounded by S q and, as above, tractions are specified on S q . Using 
conventional indicial notation, coordinates, displacements, and velocities 
in B q are d ascribed by X*, U* , an^ V*, respectively; coordinates in the de- 
formed domain B are x 1 . Large strain elasto-plastic behavior is governed 
by quasi- static rate equations of the form 
;u . pUKL- 

b " p e kl 

• I J * 

where S is the convected Kirchhoff stress rate and Ejj is the material 

derivative of Green's strain, viz.. 


* 1 if K 

p - £fv +v +U V +U V 1 
C IJ 2 lv I;J V J; I ;I K;J ;J K;J J 

IJKL 

Finally, P is a constitutive tensor which is created to coincide with 

Hooke's law for small elastic deformations but subsequently accommodates 

arbitrary work-hardening material behavior, including provision for unloading 

• , 

and reloading. The stationary principle developed by Osias [19] is then 



{P IJK 4 E 
1 IJ KL 


sK Vk V I;L WB - 


f. 


T I V J dS 


«n = o 
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where T are the traction rate;, specified on S^. With (11) one needs only to 
design an element and its interpolation functions, as above, to arrive at the 
analogue to (7) or (10). 

Osias, however, adds a c ' rther refinement which may surely be adapted 

to the ui ’nents discussed elsewhere. Because of the need to identify the 

size sc.'.le of the event*, the effective range of terms of the form (1) is also 

*This matter may be visualized as the relation between the effective radius of 
the blunted corner and that of the special element, which alters as excitation 
progresses: This radius should remain somewhat smaller than that of the element. 


considered to be a variable. Hence the approach is to define first the 
base element to be used in the host program, and then to overlay the appropriate 
counterpart of Figure 2 with special functions. For example. Figure 2 could 
be supposed to represent six-noded triangles. The sectors comprising the 
special element contain as a part* of their total interpolation functions 
displacement distributions given by 

Uj(R,0) = Rj(R/R s ,0) 0j(0-0) (no sum on I) (12) 

where R c < R is the relation between the range R c over which (12) is 
operative and the (Langrangian) radius measure of the sector; <_ 0 <_ 0£ 
is the angular coordinate, with 0 = (0^ + ^^)/2\ and the functional forms 
Rj and ^ remain to be specified. Clearly (12) is operative only for R/R g <_ 1 
and will contain (R/Rg)^; Osias observes that and its gradients must vanish 
smoothly as R/R **1. An example is suggested but not implemented in [19]; the 
point of interest here is not so much the detail as the fact that computation 
of a singularity may proceed in the context of geometric non-linearity. It 
is moreover evident that material and geometric non-linearities may be treated 
together [18,19] as a more physically realistic model than the more 
familiar linear elastic case. 

It may also be remarked that the presen-t approach admits a feature in 
modeling not elsewhere considered. We know from a variety of sources that 

plastic flow and blunting occur at a sharp, re-entrant corner, but that these 

# 

events tend to have a directional character. That is, for example, the degree 

of yield varies with 9(or 0). It is therefore not obvious that the exponent 

q in (1) or (2) should be invariant with respect to the angular coordinate. 

The present formulation allows the exponent to vary so that, for a crack, the 

relatively inactive zones tending to lie ahead of the crack and along its flanks 

*See Appendix 1 for one means for overlaying regular behavior with a special 
distribution. 
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may more realistically be modeled. Osias accounts for this situation [19], as 
did we originally [7], but there is some computational cost involved in having 
a larger number of variables in the statement (7b). Nonetheless, the analyst 
retains the option to examine directionality of non-linear behavior which, 
for some circumstances, may prove useful [15]. 

AXI SYMMETRY 

The very familiarity of linear elasticity has the virtue of facilitating 
consideration of more complex circumstances. So far the discussion has been 
in the context of planar problems owing largely to Williams's original find- 
ings. It is known that non-planar crack problems have been solved [20], and 
the question naturally arises whether the structure of this singularity re- 
lates to the planar form. While correspondence has been demonstrated for 
a crack, there is evidence *o suggest that the intensity may depend on far- 
field geometry or other features not immediate to the crack's tip or edge [21]. 

An analogue to Williams's planar eigenanalyses would be the obvious step 
for the axisymmetric case. Consideration ought not to be limited to cracks, 

t 

however, as there are other configurations of technical interest. These in- 
clude piping flanges, step changes in shaft diameter, and sudden thickness 
variations in thick shells. Unfortunately, the Williams type of analysis 

does not go forward as in the case of planar behavior, because the governing 

# 

equations do not admit a product solution [22], These equations, for the 
configuration shown in Figure 7, derive from the simple transformation 

r = R + p cos (V + , z = Z + p sin + i|> Q ) (13) 

in the (p,9,40 coordinates, displacements are (£.n,p); the strains are 


1 


I 


] 


I 
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e p * 3£/3p 

e Q = (3n/39 + Cr p + cr^/p)/r 
* Oc/34» + 0/P 

Y©^ = (3n/3<lO/p + (35/30 - nr^/p)/r 

Y^p = 3?/3p + 0£/3<l>-0/p 

y = 3n/3p + (3^/30-nr )/r 
po P 


and equilibrium is written as 

80() /3p ♦ (a p -a e )r p /r * (o p -o^)/p * Ot pe /39)/r 

* * t pt (yp)/r » 0 

3Xp 0 /3p + (1/p + 2r p /r)T pQ + (3a Q /3e)/r + (3x 0 ^/3ip3/p 

* 2T e* ( V p)/r " 0 

3T p / 3 ° * tr p /r * 2/o)T P * * (3T e* /39)/r * t3 V 3H/p 

* ( V°e )C V p)/r ’ 0 

« 

In (14) and (15) we use the notation 

r p = cos (i|>+i|/ o ), r^/p = - sin (♦ + 4» Q ) 

For axisymmetry, 3/30 is a null operator and (14) (15) decouple to two sets, 

one for extension and one for torsion in which 

0 

£ = €(p, 40, n = 0, 5 = 5 (p,i| 0 and K • 0, n = n(p,t). 5 = 0 

For extension, functions of the Neuber-Papkovich type may be introduced 
after consideration of the coordinate transformation [22]; denoting these 


as a) and d, we find 


X I [ I I ! ) 1 
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2 u5 * -3w/3p + (3-4v)tf - [p+R cos(ij/+iJi o ) ]3$/3p 

2y; = -(l/p)3w/3^ - [l + (R/p)cos(ii»+ii/ o )]3«>/94» + [(R/p)sinO+(J> 0 ) ]$ 

The corresponding differential equations are decoupled but awkward so that 
we introduce 
A = (/r)a) 

0 * (/r)«j 
to find 

3 2 A/3p 2 + (1/p) 3A/3p + (l/p 2 )3 2 A/3^ 2 ♦ (l/4r 2 )A = 0 

3 2 0/3p 2 + (1/p) 30/3p + (1 /p 2 )3 2 0/3i^ 2 (16) 

♦ [(l/4-r p 2 )/r 2 -l/p 2 ]0 = 0 

as the equations of interest. (Forms similar to the first of (16) obtain 
for torsion.) The presence of both p and r in the denominators of the 
coefficient in the various terms of (14), (IS), and (16) precludes the "sep- 
aration of variables" approach inherent in (8) or (9), and the template for 
solution provided by Williams fails*. It is, of course, possible to assume 
p « R and integrate (16) in terms of Bessel functions , seeking asymptotic 
forms for small p. Such an approach, however, is objectionable for two 
reasons. The solution is approximate and equilibrium is not fully satisfied. 
Further, the degree of approximation remains unestablished and there may be 
errors at the boundary itself. Second, such a solution is not useful in 
procedures already developed for such problems including ordinary special 
elements, e.g., [2,3], edge functions [4], and collocation [5]; one is left 
with finite elements as in [21]. 

*We have also examined formulations including the Galerkin vector and its 
special form. Love’s strain function; Southwell functions, both as originally 
stated [23] and as modified by Zak [24]; and the Maxwell-Morera functions. 

Of these, the procedure shown is the most promising. 


If we must settle for approximate solutions, it seems sensible to 
adjust the approximation to the particular problem. Thus we prefer an 
approach which is patterned after that outlined in this paper. The analogue 
to (2) becomes simply 

u = u + A.r + C.z + (E+Fil/)p^ cos (iL+ip ) - (G+H40p^ sin (4< + ^ ) 
r o 1 1 o o (1?a) 

w z = w q + D^r + B^z + (E+FiJOp^ sin + (G+HijOp^ cos (iJj+4^) 

and 

5 - u Q cos + w q sin 0 +iJj o ) + yCAj+B^p + cos [2 ] 

+ ^-(C.+D )p sin [2 (i|i+ 4> )] + (E+DiJOp^ 

1 1 1 0 (17b) 

n = -u sin (t+^ 0 ) + w q cos (<Ji+i[> 0 ) - jCC^-D^p + y( c 1 +D 1 )P cos [2 (i|H-ip Q ) ] 

- ^‘(A^-Bj) p sin [2(4 j+^ o ) ] + (G+Hip)p q 

Using (14), the strains are found and the procedure of computing and mini- 
mizing potential energy is followed, directly in analogue to the foregoing 
discussion. 

As has been observed, by many authors, probably first by Irwin [25], 
the singularity for a crack geometry is expected to be the familiar value 
q = 1/2. What is not known, however, is the result for geometries associated 
with a flange, a change of section, layered materials, and so on. That is, 
what is q(R,a,4i o )? To the extent that such information is of technical 
interest, either for use in other computations (e.g., [3,4]) or for specific 
applications, the foregoing development appears to be the first direct means 
available for establishing the structure of an elastic singularity at a 
rfe-entrant corner in an axisymmetric geometry. 

Moreover, it is evident from the work now in progress (which deals with 
non-linearities in planar situations) that the procedure carries over directly 


to non-linearities in the axisymmetric case. In the same sense that finite 
elements are used in both planar and axisymmetric analyses, extension of the 
formulation in (17) is appropriate to elasto-plastic flow [14] and to large 
deformations [18,19], For these situations, the nature of the singularity 
may be established numerically and perhaps in no other manner. 

THREE DIMENSIONS 

The numerical approach described here is fully pertinent to three- 
dimensional problems. The context is a vertex formed by intersection of 
three surfaces, and two simple examples are sketched in Figure 8. Further 
examples may be drawn from the literature in crack problems, where a crack 
intersects a free surface. The crack front may be straight or curved, and 
the intersection need not be normal to the free surface. Some possibilities 
were outlined earlier [21], and Dr. B. K. Neale has recently reported a 
series of such observations [26], From experimental sources, the inference 
of a singularity peculiar to the vertex itself is not uncommon, and Prof- 
essor E. S. Folias has made an initial analytic extraction of such behavior [27]. 
The highly resolved numerical results developed by Dr. T. A. Cruse are also 
supportitve of the notion of a vertex singularity [28]. 

Such information is appropriate to the computational approach outlined 
here. It is easy to envisage the shapes sketched in Figure 8 filled with 
polygonal cones whose apices all coincide at the vertex of the corner. A 
radial coordinate from that point becomes p in (1). The question then becomes 
how the displacements vary with this radial coordinate. There will of course 
be rigid motion, and one must also allow for constant strains in the usual 
manner. Obviously, however, behavior of the sort suggested by (1) can intrude 
as well; information of the sort outlined above suggests that inclusion of 
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terms of the form (1) is most appropriate. As a result, one may infer a 
process of modeling which will explicitly account for a vertex singularity. 
Furthermore, one may study behavior along an edge (intersection of two sur- 
faces) by the same means and examine, among other points, structure of the 
singularity as a function of the shape of the edge. 

It would be tedious here to write the equations of elasticity in spher- 
ical coordinates centered at the vertex, and to expand them for elementary 
interpolation functions for, say, a tetrahedron or a pyramid*. Obviously 
enough, such development is straightforward if lengthy to carry out, and it 
could be implemented in terms of code. We prefer instead to emphasize that 
the discussion following (9) is immediately applicable. Minimization of 
the potential energy yields an approximate statement of equilibrium (sector- 
wise as opposed to point-wise) and effectively exact boundary conditions. 

To the extent that vhe exponent on the radial coordinate is an artifact of 
these boundary conditions and not the field equations, the analyst may 
anticipate determination of the vertex singularity to a considerable degree 
of accuracy. Moreover, the singularity computation is direct and explicit, 
and it should reveal the structure of the singularity without need for 
extensive data reduction. 

CONCLUDING REMARKS 

# - ? 

We have sought to describe an approach for singularity computations 
which is based on a primary concept in mechanics. The procedure replicates 
the conditions used for the same purpose in more formal analysis, as applied 
to planar configurations, but is in no manner limited to linear or elastic 
or small strain or two-dimensional situations. Rather, the range of cir- 


*That is to say, a triangular cone or a rectangular cone. 


cumstances for which the approach may be used appears to be unlimited, and 
it evidently carries over to other classes of problems such as heat flow. 
As a result, there is now an opportunity to attack problems which so far 
have proven difficult in the sense that details of the singularity’s 
structure are not yet well articulated. While this approach requires more 
extensive preparation prior to implementation and additional computation 
costs are encountered, the directness of obtaining results needed for 
certain research objectives would appear to make this approach worthy of 
use in serious singularity computations. 
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APPENDIX - ADDITIONAL INTERPOLATION FUNCTIONS 


The element formulation in (2) exemplifies the basic approach needed 
to exploit (1) in the sense of an undetermined exponent. By no means, how- 
ever, is the procedure limited to the primitive form in (2). Much more sub- 
stance can be incorporated in the interpolation functions, but there are 
evident constraints to be observed. Noting that (1) derives from a certain 
physical sense of the solution's behavior, one must guard against the use of 
coordinates whose physical interpretation is not clearly understood. Thus, 
for example, use of isoparametric elements is problematic unless the internal 
coordinate(s) can be directly identified with the radial coordinate of interest. 
Indeed, Freese and Tracey [29] recently noted the sensitivity of elastic be- 
havior in crack problems to overall element shape. 

It appears therefore advisable to work with physical coordinates when 
formulating an element, or to use an isoparametric form in which the physical 
coordinates are manifest. Marino, in work soon to be documented, employs 
the latter method with evident success [15]. An example of the former is 
suggested in Figure 9 and a more complex version of (2) as follows: 


2 2 

u * u q + A^x + A^x + C^y + + Exy 

+ [G 0 + G 1 (e-?) + G 2 (e-F) 2 ]p p cos e 

0 

- [H q + H^e-9) + H 2 (6-F) 2 ]p q sin 6 
2 • 2 

v * v q ♦ D^x + D 2 x + Bjy + B 2 y + Fxy 

+ [G o + G 1 (0-0) + G 2 (6-F) 2 ]p P sin 0 
+ [H + H.( 0 - 0 ) + a,(0-e) 2 ]p q cos 0 

Ol L 

where 0 = (0^+0 2 )/2. It follows that 


(18a) 
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u ■ u cos 0 ♦ v sin 0 + r(A.+i> )p ♦ ^(A.-B }p cos 20 
p o o 2 ) 1 2 1 1 

♦ k V D l)P sin 2i ♦ j(A 2 +C 2 +F)p 2 cos 0 

+ y(A 2 -C 2 -F)p 2 cos 0 cos 20 ♦ j(B 2 +D 2 *E) p 2 sin 0 

+ |(B 0 -D,-E)p 2 sin 0 cos 20 ♦ [G • G ( 0 -?) ♦ G,(0-0) 2 ]p P 

2 2 2 oi i (18b) 

v. * -u sin 0 + v cos 0 - 4(C.-D,)p + icC.+D^p cos 20 
0o o 2 1 1' 2 1 V 

- yCAj-B^p sin 20 - ±(A 2 +C 2 -F)p 2 sin 0 

- y(A 2 -C 2 -F)p 2 sin 0 cos 20 + y(B 2 +D 2 -E)p 2 cos 0 

- ^(B 2 -D 2 -E)p 2 cos Q cos 20 ♦ [H q +H (0-6) ♦ H 2 (9-e) 2 ]p q 

Note that two exponents are allowed in this representation, one for each of 
the polar components of displacement. The remainder of the expression is 
based on the notion of a linear strain element (LSE). In that context, note 
that an eight-noded element obtains by omitting the G 2 and H 2 terms in (18b) 
or perhaps the E and F terms in (18a). We could, in addition, let th. expo- 
nents vary in the angular direction, i.e., p,q « p(9),q(0), but care must 
be exercised to avoid intersector incompatabilities . 

A feature of re-entrant corner and crack problems, sometimes overlooked, 
is that circumferential gradients of strain and stress can be greater than 
radial gradients. This may be observed most simply by plotting contours of 
selected field quantities, as in [30], and noting behavior in the respective 
directions. One may also use this as a basis for sizing the special element, 
or its effective radial range. Alternatively, an energy basis may be used 
as described in the main part <i the paper. 

Other formulations may be generated in this intuitive manner, each of 
which should be evaluated for its efficacy with respect to the problem type 
anticipated. It is not so much our purpose here to devise an optimal element - 
that ) Ution appears to be problem dependent - but we do note that formulation 
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may proceed independently of the overall conceptual framework. In that con- 
text, we should note the procedure termed ’’singularity programming" devised 
by Emery and his colleagues {31], This formulation involves superposing 
special or irregular behavior on an established element by writing 

(u> - [R](A) ♦ k{S> (19) 

Here, (u) has the same meaning as in (3); [R] is a matrix of regular, typ- 
ically analytic, interpolation functions as in the CSE, LSE, or isopara- 
metric representation; (A) is a vector of constants; k is a scale factor; 
and {S} is a vector of interpolation functions which contain the exponent q 
and represent the special aspects of the element. Evaluating (u) at the 
nodes leads to the vector (uh (19) is solved for (A) and rewritten as 
(u) - (N](u) ♦ k{S) 

where 

(NJ . [R] [R]" 1 
{S> * <S> - IN] (SJ 

Clearly (S) vanishes at the nc ?r; it does, however, interact with the nodal 
excitation as may be seen- by carrying through development of the stiffness 
equation. The formulation is interesting in that, vMle it allows determina- 
tion of the exponent (s) in { S} , the magnitude of the special behavior is 
scaled by the single factor k. Hence one must either be more specific as 
to the angular variation, relative to what is permitted in (18), or settle 
for possibly a stiffer element. In addition, since the potential energy is 
minimized with respect to k, there are entries in the stiffness equation 
beyond the usual terms, whereas the formulations (2) and (18) do not present 
such an inconsistency. Whatever admissible formulation is used (see also 
[32]), the option of adjusting the singular structure remains available and 
the analyst must choose ultimately in terms of the problem class he faces. 
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y, v 0 < p < Zp t 



Figure 1, showing a special element comprising twenty sectors and 
a typical sector with details of coordinates, displacement components, 

and dimensions. 
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Figure 2, showing a special element embedded in regular elements, 
CSE formulation. 
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Figure 3, showing radial variation of octahedral stress 
normalized on limit value (at which yield initiates) at 
0 = 90 deg for four load levels; data from problem reported 
in [11]. 


(u-u 0 )/A x 10 2 at 9 =90 deg 
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Figure 5, showing angular variation of octahedral stress 
(normalized on limit value) at r/a ^ 0.008 for four load 
levels; data from problem reported in [11] (note 4x magni- 
fication of elastic results). 


(u-u 0 )/A xIO at r/a~0.008 
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Figure 6, showing angular variation of normalized nodal displacements at 
r/a ^ 0.008 for four load levels; from [11]. 
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Figure 7, showing coordinates p , centered at notch vertex in 
an axisymmetric coordinate system r,z. Note that origin for p 
is at R, Z and that i|j = 0 bisects (interior) notch angle 2a. 



Figure 8, showing cusp-like and re-entrant vertices in three 
dimensions, formed by intersection of three smooth surfaces. 
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Figure 9, showing an example of a sector in a more refined 
special element. Details are presumed to follow the pattern 
shown in Figure 1. 


